# score calculation & merge threshold calculation
# 2014/01/29

treatmentDir <- commandArgs(trailingOnly =T)[1]
controlDir <- commandArgs(trailingOnly =T)[2]
ratio <- as.numeric(commandArgs(trailingOnly =T)[3])
FileScore <- commandArgs(trailingOnly =T)[4]
FileMergeThreshold <- commandArgs(trailingOnly =T)[5]
FileAverageTreatment <- commandArgs(trailingOnly =T)[6]

### title
cat("#chr", "pos", "treatment", "control", "substract_nor", "logPval", "fold_change\n", sep="\t", file=FileScore, append=FALSE)



readSGR = function(chr){
	file_tr = paste(treatmentDir, "/", chr, ".sgr", sep="")
	data_tr = read.table(file_tr, header=F, sep="\t")
	tr = data_tr[,2]
	file_co = paste(controlDir, "/", chr, ".sgr", sep="")
	data_co = read.table(file_co, header=F, sep="\t")
	co = data_co[,2] * ratio	

	maxIndex = length(tr)
	if(length(co) < maxIndex){
		maxIndex = length(co)
	}

	substract = tr[1:maxIndex] - co[1:maxIndex]
	fold_change = (tr[1:maxIndex] + 1) / (co[1:maxIndex] + 1)
	
	fr <- data.frame(chr=rep(chr,maxIndex), pos=data_tr[1:maxIndex,1], treatment=tr[1:maxIndex], control=co[1:maxIndex], substract=substract, fold_change=fold_change)
	
	# return frame
	fr
}

entireData <- readSGR("I")
entireData <- rbind(entireData, readSGR("II"))
LengthChrIChrII = length(entireData$substract)
entireData <- rbind(entireData, readSGR("III"))
entireData <- rbind(entireData, readSGR("MTR"))
entireData <- rbind(entireData, readSGR("MT"))

attach(entireData)

### normalize using without chrIII
IndexChrIChrII = 1:LengthChrIChrII
substract_norm = substract / mean(treatment[IndexChrIChrII])

cat(mean(treatment[IndexChrIChrII]), file=FileAverageTreatment)

logPval = pnorm(substract, mean(substract[IndexChrIChrII]), sd(substract[IndexChrIChrII]), lower.tail=FALSE, log.p = TRUE)
output <- data.frame(chr=chr, pos=pos, treatment=treatment, control=round(control,3), substract_norm=round(substract_norm, 5), logPval=round(logPval,5), fold_change=round(fold_change,5))
write.table(output, file=FileScore, append=TRUE, quote=F, sep="\t", row.names=F, col.names=F)

threshold = qpois(0.001, mean(treatment[IndexChrIChrII], na.rm=TRUE), lower.tail=FALSE)
#message <- paste("\nThreshold for merge is ", as.character(threshold), "\n", sep="")
#cat(message)

cat(threshold, file=FileMergeThreshold)
